In [1]:
import numpy as np
import plotly.graph_objects as go
from plotly.subplots import make_subplots
from scipy import constants

a0 = constants.physical_constants['Bohr radius'][0]
Angstrom = constants.angstrom

Here is a function to calculate the spherical coordinates r , ϕ , θ of a point with Cartesian coordinates x , y , z .

In [2]:
def r_phi_theta(x, y, z):
    r = np.sqrt(x**2 + y**2 + z**2)
    phi = np.arctan2(y, x) + np.pi
    theta = np.arccos(z/r)
    return r, phi, theta

The next cell uses numpy.mgrid to create a 3d grid x y z grid spanning [ − 20 , 20 ] Angstroms in each dimension. The grid has 48

3¶

110592 gridpoints. More gridpoints means better-looking graphics but also longer compute times. mgrid gives us three individual arrays: one contains the x coordinates of all points in our 3d grid, another has the y coordinates, and the third has the z coordinates.

The second line of code uses the r_phi_theta function defined above to calculate three arrays that give us the same 3d grid as above, but in spherical coordinates.

You will notice how the output values of both mgrid and r_phi_theta() are "unpacked" and assigned to three separate values. Non-unpacked syntax such as

grid = r_phi_theta(X, Y, Z) would return a tuple containing the three spherical coordinate grids.

In [3]:
X, Y, Z = np.mgrid[-20*a0:20*a0:48j, -20*a0:20*a0:48j, -20*a0:20*a0:48j]
r_grid, phi_grid, theta_grid = r_phi_theta(X, Y, Z)

The next cell contains a probability function that performs the operation r 2 ψ ∗ ( r , ϕ , θ ) ψ ( r , ϕ , θ ) . The second line normalizes the probability function in a way that makes it easy to plot. The requirement for the probability function is that ∫ ∞

r¶

0 ∫ 2 π

ϕ¶

0 ∫ π

θ¶

0 r 2 ψ ( r , ϕ , θ ) ∗ ψ ( r , ϕ , θ ) d r d ϕ d

θ¶

1. However, since the volume of the hydrogen atom is so tiny (recall that the Bohr radius is less than an Angstrom), the probability densities can be huge, of order 10 26 m − 3 . We will put our probability densities on a scale of 0 to 1 just so the numbers on our colorbars won't look insane.

In [4]:
def probfunc(r, wavefunc):
    p = r**2 * np.conj(wavefunc) * wavefunc
    return p.real / np.max(p.real)
# The .real syntax picks off the real part of a complex number. Even though the 
# imaginary part of the probability is zero, variable p is still represented in
# the computer as a complex number. Our plotting routines will only work on
# arrays of real numbers.

The next cell defines two wave functions: n1l0m0(r) for n , l ,

m¶

1 , 0 , 0 and n2l0m0(r) for n , l ,

m¶

2 , 0 , 0 .

Question to answer:¶

How come these wave functions do not depend on ϕ and θ ?

Since the atom we are dealing with is spherical, it has symmetry at all angles of ϕ and θ thus determining the electron potential energy only depends on r, the distance from the nucleus. The symmetry removes any dependence on ϕ and θ.¶

In [ ]:
n1l0m0 = lambda r: np.exp(-r/a0) / (np.sqrt(np.pi) * a0**1.5)

n2l0m0 = lambda r: (2 - r/a0) * np.exp(-r/(2*a0)) / (4 * np.sqrt(2*np.pi) * a0**1.5)

The next cell finds the probability functions associated with n , l ,

m¶

1 , 0 , 0 and n , l ,

m¶

2 , 0 , 0 :

In [ ]:
prob_n1l0m0 = probfunc(r_grid, n1l0m0(r_grid))
prob_n2l0m0 = probfunc(r_grid, n2l0m0(r_grid))

Now we use plotly to plot isosurfaces of our probability functions An isosurface is like a contour line on a map, except in 3d: it's a surface over which the value of the dependent variable is constant. These visualizations will each have eight semi-transparent isosurfaces. Darker surfaces correspond to higher probability densities. You can click and drag to rotate, enlarge, or shrink your visualization.

You will see from the graphics that n , l ,

m¶

1 , 0 , 0 and n , l ,

m¶

2 , 0 , 0 give spherically symmetric probability functions, as expected.

n =1 ,

ℓ¶

0 , m

l¶

0

In [ ]:
# n = 1, l = 0, m = 0
fig = go.Figure(data=go.Isosurface(
    x=X.flatten() / Angstrom, # x values of the grid in Angstroms
    y=Y.flatten() / Angstrom, # y values of the grid in Angstroms
    z=Z.flatten() / Angstrom, # z values of the grid in Angstroms
    value=prob_n1l0m0.flatten(), # independent variable
    isomin=0.05, # Minimum normalized probability density to render in an isosurface
    isomax=0.95, # Maximum normalized probability density to render in an isosurface
    opacity=0.4, # Set a low opacity so each surface is partially transparent
    colorscale='Plotly3_r', # Nice-looking color table
    surface_count=8, # number of isosurfaces to plot (2 by default: only min and max)
    colorbar_nticks=8, # colorbar ticks correspond to isosurface values
    caps=dict(x_show=False, y_show=False, z_show=False))
               )

# Change axis lables and make the plot larger than the default
fig.update_layout(scene = dict(
                    xaxis_title='x (Angstroms)',
                    yaxis_title='y (Angstroms)',
                    zaxis_title='z (Angstroms)'),
                    width=700,
                    margin=dict(r=10, b=10, l=10, t=10))

fig.show()
In [ ]:
# n = 2, l = 0, m = 0
fig = go.Figure(data=go.Isosurface(
    x=X.flatten() / Angstrom, # x values of the grid in Angstroms
    y=Y.flatten() / Angstrom, # y values of the grid in Angstroms
    z=Z.flatten() / Angstrom, # z values of the grid in Angstroms
    value=prob_n2l0m0.flatten(), # independent variable
    isomin=0.05, # Minimum normalized probability density to render in an isosurface
    isomax=0.95, # Maximum normalized probability density to render in an isosurface
    opacity=0.4, # Set a low opacity so each surface is partially transparent
    colorscale='Plotly3_r', # Nice-looking color table
    surface_count=8, # number of isosurfaces to plot (2 by default: only min and max)
    colorbar_nticks=8, # colorbar ticks correspond to isosurface values
    caps=dict(x_show=False, y_show=False, z_show=False)
    ))

# Change axis lables and make the plot larger than the default
fig.update_layout(scene = dict(
                    xaxis_title='x (Angstroms)',
                    yaxis_title='y (Angstroms)',
                    zaxis_title='z (Angstroms)'),
                    width=700,
                    margin=dict(r=10, b=10, l=10, t=10))

fig.show()

Question: which of the two probability functions plotted above is more centrally condensed, that is, a higher probability of the electron being close to the proton?¶

Where n=1, the first plot has a higher probability of the electron being close to the proton because it has less points that extend out to further angstroms¶

Slices¶

Our isosurface renderings might not help us visualize a probability density that is not monotonically increasing as a function of radius: a low probability region (light pink) might be concealed insided a higher-probability region (purple). To visualize the radial structure of our probability functions in more detail, we add a slice to each isosurface rendering. The slice I chose is parallel to the y − z plane and centered at

x¶

0 . See if you can spot the one-line change to the code that adds the slice.

n =1 ,

ℓ¶

0 ,

m¶

0

In [ ]:
# n = 1, l = 0, m = 0
fig = go.Figure(data=go.Isosurface(
    x=X.flatten() / Angstrom, # x values of the grid in Angstroms
    y=Y.flatten() / Angstrom, # y values of the grid in Angstroms
    z=Z.flatten() / Angstrom, # z values of the grid in Angstroms
    value=prob_n1l0m0.flatten(), # independent variable
    isomin=0.05, # Minimum normalized probability density to render in an isosurface
    isomax=0.95, # Maximum normalized probability density to render in an isosurface
    opacity=0.4, # Set a low opacity so each surface is partially transparent
    colorscale='Plotly3_r', # Nice-looking color table
    surface_count=8, # number of isosurfaces to plot (2 by default: only min and max)
    colorbar_nticks=8, # colorbar ticks correspond to isosurface values
    slices_x=dict(show=True, locations=[0]),
    caps=dict(x_show=False, y_show=False, z_show=False))
               )

# Change axis lables and make the plot larger than the default
fig.update_layout(scene = dict(
                    xaxis_title='x (Angstroms)',
                    yaxis_title='y (Angstroms)',
                    zaxis_title='z (Angstroms)'),
                    width=700,
                    margin=dict(r=10, b=10, l=10, t=10))

fig.show()

Does the slice shown above match your expectations of the probability function?

The slice does match the expectations of the probability function because it is mostly centrally condensed.¶

n =2 ,

ℓ¶

0 ,

m¶

0

In [ ]:
# n = 2, l = 0, m = 0
fig = go.Figure(data=go.Isosurface(
    x=X.flatten() / Angstrom, # x values of the grid in Angstroms
    y=Y.flatten() / Angstrom, # y values of the grid in Angstroms
    z=Z.flatten() / Angstrom, # z values of the grid in Angstroms
    value=prob_n2l0m0.flatten(), # independent variable
    isomin=0.05, # Minimum normalized probability density to render in an isosurface
    isomax=0.95, # Maximum normalized probability density to render in an isosurface
    opacity=0.4, # Set a low opacity so each surface is partially transparent
    colorscale='Plotly3_r', # Nice-looking color table
    surface_count=8, # number of isosurfaces to plot (2 by default: only min and max)
    colorbar_nticks=8, # colorbar ticks correspond to isosurface values
    caps=dict(x_show=False, y_show=False, z_show=False),
    slices_x=dict(show=True, locations=[0]),
    ))

# Change axis lables and make the plot larger than the default
fig.update_layout(scene = dict(
                    xaxis_title='x (Angstroms)',
                    yaxis_title='y (Angstroms)',
                    zaxis_title='z (Angstroms)'),
                    width=700,
                    margin=dict(r=10, b=10, l=10, t=10))

fig.show()

Based on what you see in the slices, how are the probability functions for

n¶

1 ,

ℓ¶

0 ,

m¶

0 and

n¶

2 ,

ℓ¶

0 ,

m¶

0 different?

The center of the probabilty function is less dense with higher probabilities in the n=2 function. The center has rings of different probabilities which make is less likely to be very close to the nucleus until it gets to a certain range r away.¶

Try

n¶

2 ,

ℓ¶

1 ,

m¶

0

In [ ]:
n2l1m0 = lambda r, theta: (r/a0) * np.exp(-r/(2*a0)) * np.cos(theta) \
                          / (4 * np.sqrt(2*np.pi) * a0**1.5)
In [ ]:
# n = 2, l = 1, m = 0
prob_n2l1m0 = probfunc(r_grid, n2l1m0(r_grid, theta_grid))

fig = go.Figure(data=go.Isosurface(
    x=X.flatten() / Angstrom, # x values of the grid in Angstroms
    y=Y.flatten() / Angstrom, # y values of the grid in Angstroms
    z=Z.flatten() / Angstrom, # z values of the grid in Angstroms
    value=prob_n2l1m0.flatten(), # independent variable
    isomin=0.05, # Minimum normalized probability density to render in an isosurface
    isomax=0.95, # Maximum normalized probability density to render in an isosurface
    opacity=0.4, # Set a low opacity so each surface is partially transparent
    colorscale='Plotly3_r', # Nice-looking color table
    surface_count=8, # number of isosurfaces to plot (2 by default: only min and max)
    colorbar_nticks=8, # colorbar ticks correspond to isosurface values
    caps=dict(x_show=False, y_show=False, z_show=False),
    # slices_x=dict(show=True, locations=[0]),
    ))

# Change axis lables and make the plot larger than the default
fig.update_layout(scene = dict(
                    xaxis_title='x (Angstroms)',
                    yaxis_title='y (Angstroms)',
                    zaxis_title='z (Angstroms)'),
                    width=700,
                    margin=dict(r=10, b=10, l=10, t=10))

fig.show()

n =2 ,

ℓ¶

1 ,

m¶

1

In [ ]:
n2l1m1 = lambda r, phi, theta: (r/a0) * np.exp(-r/(2*a0)) * np.sin(theta) \
                               * np.exp(1j * phi) / (8 * np.sqrt(3*np.pi) * a0**1.5)
In [ ]:
# n = 2, l = 1, m = 1
prob_n2l1m1 = probfunc(r_grid, n2l1m1(r_grid, phi_grid, theta_grid))

fig = go.Figure(data=go.Isosurface(
    x=X.flatten() / Angstrom, # x values of the grid in Angstroms
    y=Y.flatten() / Angstrom, # y values of the grid in Angstroms
    z=Z.flatten() / Angstrom, # z values of the grid in Angstroms
    value=prob_n2l1m1.flatten(), # independent variable
    isomin=0.05, # Minimum normalized probability density to render in an isosurface
    isomax=0.95, # Maximum normalized probability density to render in an isosurface
    opacity=0.4, # Set a low opacity so each surface is partially transparent
    colorscale='Plotly3_r', # Nice-looking color table
    surface_count=8, # number of isosurfaces to plot (2 by default: only min and max)
    colorbar_nticks=8, # colorbar ticks correspond to isosurface values
    caps=dict(x_show=False, y_show=False, z_show=False),
    # slices_x=dict(show=True, locations=[0]),
    ))

# Change axis lables and make the plot larger than the default
fig.update_layout(scene = dict(
                    xaxis_title='x (Angstroms)',
                    yaxis_title='y (Angstroms)',
                    zaxis_title='z (Angstroms)'),
                    width=700,
                    margin=dict(r=10, b=10, l=10, t=10))

fig.show()
In [ ]:
# n = 2, l = 1, m = 1
fig = go.Figure(data=go.Isosurface(
    x=X.flatten() / Angstrom, # x values of the grid in Angstroms
    y=Y.flatten() / Angstrom, # y values of the grid in Angstroms
    z=Z.flatten() / Angstrom, # z values of the grid in Angstroms
    value=prob_n2l1m1.flatten(), # independent variable
    isomin=0.05, # Minimum normalized probability density to render in an isosurface
    isomax=0.95, # Maximum normalized probability density to render in an isosurface
    opacity=0.4, # Set a low opacity so each surface is partially transparent
    colorscale='Plotly3_r', # Nice-looking color table
    surface_count=8, # number of isosurfaces to plot (2 by default: only min and max)
    colorbar_nticks=8, # colorbar ticks correspond to isosurface values
    caps=dict(x_show=False, y_show=False, z_show=False),
    slices_x=dict(show=True, locations=[0]),
    ))

# Change axis lables and make the plot larger than the default
fig.update_layout(scene = dict(
                    xaxis_title='x (Angstroms)',
                    yaxis_title='y (Angstroms)',
                    zaxis_title='z (Angstroms)'),
                    width=700,
                    margin=dict(r=10, b=10, l=10, t=10))

fig.show()
In [ ]:
# n = 2, l = 1, m = 1
fig = go.Figure(data=go.Isosurface(
    x=X.flatten() / Angstrom, # x values of the grid in Angstroms
    y=Y.flatten() / Angstrom, # y values of the grid in Angstroms
    z=Z.flatten() / Angstrom, # z values of the grid in Angstroms
    value=prob_n2l1m1.flatten(), # independent variable
    isomin=0.05, # Minimum normalized probability density to render in an isosurface
    isomax=0.95, # Maximum normalized probability density to render in an isosurface
    opacity=0.4, # Set a low opacity so each surface is partially transparent
    colorscale='Plotly3_r', # Nice-looking color table
    surface_count=8, # number of isosurfaces to plot (2 by default: only min and max)
    colorbar_nticks=8, # colorbar ticks correspond to isosurface values
    caps=dict(x_show=False, y_show=False, z_show=False),
    slices_z=dict(show=True, locations=[0]),
    ))

# Change axis lables and make the plot larger than the default
fig.update_layout(scene = dict(
                    xaxis_title='x (Angstroms)',
                    yaxis_title='y (Angstroms)',
                    zaxis_title='z (Angstroms)'),
                    width=700,
                    margin=dict(r=10, b=10, l=10, t=10))

fig.show()

A fancier wave function:

n¶

3 ,

ℓ¶

2 ,

m¶

1

In [ ]:
n3l2m1 = lambda r, phi, theta: (r**2/a0**2) * np.exp(-r/(3*a0)) * np.sin(theta) * \
                               np.cos(theta) * (np.cos(phi) + 1j * np.sin(phi)) / \
                               (81 * np.sqrt(np.pi) * a0**1.5)
In [ ]:
# n = 3, l = 2, m = 1
prob_n3l2m1 = probfunc(r_grid, n3l2m1(r_grid, phi_grid, theta_grid))

fig = go.Figure(data=go.Isosurface(
    x=X.flatten() / Angstrom, # x values of the grid in Angstroms
    y=Y.flatten() / Angstrom, # y values of the grid in Angstroms
    z=Z.flatten() / Angstrom, # z values of the grid in Angstroms
    value=prob_n3l2m1.flatten(), # independent variable
    isomin=0.05, # Minimum normalized probability density to render in an isosurface
    isomax=0.95, # Maximum normalized probability density to render in an isosurface
    opacity=0.4, # Set a low opacity so each surface is partially transparent
    colorscale='Plotly3_r', # Nice-looking color table
    surface_count=8, # number of isosurfaces to plot (2 by default: only min and max)
    colorbar_nticks=8, # colorbar ticks correspond to isosurface values
    caps=dict(x_show=False, y_show=False, z_show=False),
    ))

# Change axis lables and make the plot larger than the default
fig.update_layout(scene = dict(
                    xaxis_title='x (Angstroms)',
                    yaxis_title='y (Angstroms)',
                    zaxis_title='z (Angstroms)'),
                    width=700,
                    margin=dict(r=10, b=10, l=10, t=10))

fig.show()
In [ ]:
# n = 3, l = 2, m = 1
prob_n3l2m1 = probfunc(r_grid, n3l2m1(r_grid, phi_grid, theta_grid))

fig = go.Figure(data=go.Isosurface(
    x=X.flatten() / Angstrom, # x values of the grid in Angstroms
    y=Y.flatten() / Angstrom, # y values of the grid in Angstroms
    z=Z.flatten() / Angstrom, # z values of the grid in Angstroms
    value=prob_n3l2m1.flatten(), # independent variable
    isomin=0.05, # Minimum normalized probability density to render in an isosurface
    isomax=0.95, # Maximum normalized probability density to render in an isosurface
    opacity=0.4, # Set a low opacity so each surface is partially transparent
    colorscale='Plotly3_r', # Nice-looking color table
    surface_count=8, # number of isosurfaces to plot (2 by default: only min and max)
    colorbar_nticks=8, # colorbar ticks correspond to isosurface values
    caps=dict(x_show=False, y_show=False, z_show=False),
    slices_x=dict(show=True, locations=[0]),
    ))

# Change axis lables and make the plot larger than the default
fig.update_layout(scene = dict(
                    xaxis_title='x (Angstroms)',
                    yaxis_title='y (Angstroms)',
                    zaxis_title='z (Angstroms)'),
                    width=700,
                    margin=dict(r=10, b=10, l=10, t=10))

fig.show()

Your assignment: make a 3d isosurface rendering of a probabilty function that is not plotted above.¶

You may use one of the wave functions tabulated at the top of this assignment or find a more complex one if you're interested in intricate shapes.

n =3 , ℓ =0 , m =0

In [5]:
n3l0m0 = lambda r: np.exp(-r/(2*a0)) * (27 - (18*r/a0) + (2*r**2/a0**2)) / (81 * np.sqrt(3 * np.pi) * a0**1.5)
In [6]:
# n = 3, l = 0, m = 0
prob_n3l0m0 = probfunc(r_grid, n3l0m0(r_grid))

fig = go.Figure(data=go.Isosurface(
    x=X.flatten() / Angstrom, # x values of the grid in Angstroms
    y=Y.flatten() / Angstrom, # y values of the grid in Angstroms
    z=Z.flatten() / Angstrom, # z values of the grid in Angstroms
    value=prob_n3l0m0.flatten(), # independent variable
    isomin=0.05, # Minimum normalized probability density to render in an isosurface
    isomax=0.95, # Maximum normalized probability density to render in an isosurface
    opacity=0.4, # Set a low opacity so each surface is partially transparent
    colorscale='Plotly3_r', # Nice-looking color table
    surface_count=8, # number of isosurfaces to plot (2 by default: only min and max)
    colorbar_nticks=8, # colorbar ticks correspond to isosurface values
    caps=dict(x_show=False, y_show=False, z_show=False),
    ))

# Change axis lables and make the plot larger than the default
fig.update_layout(scene = dict(
                    xaxis_title='x (Angstroms)',
                    yaxis_title='y (Angstroms)',
                    zaxis_title='z (Angstroms)'),
                    width=700,
                    margin=dict(r=10, b=10, l=10, t=10))

fig.show()
In [7]:
# n = 3, l = 0, m = 0
prob_n3l2m1 = probfunc(r_grid, n3l0m0(r_grid))

fig = go.Figure(data=go.Isosurface(
    x=X.flatten() / Angstrom, # x values of the grid in Angstroms
    y=Y.flatten() / Angstrom, # y values of the grid in Angstroms
    z=Z.flatten() / Angstrom, # z values of the grid in Angstroms
    value=prob_n3l0m0.flatten(), # independent variable
    isomin=0.05, # Minimum normalized probability density to render in an isosurface
    isomax=0.95, # Maximum normalized probability density to render in an isosurface
    opacity=0.4, # Set a low opacity so each surface is partially transparent
    colorscale='Plotly3_r', # Nice-looking color table
    surface_count=8, # number of isosurfaces to plot (2 by default: only min and max)
    colorbar_nticks=8, # colorbar ticks correspond to isosurface values
    caps=dict(x_show=False, y_show=False, z_show=False),
    slices_x=dict(show=True, locations=[0]),
    ))

# Change axis lables and make the plot larger than the default
fig.update_layout(scene = dict(
                    xaxis_title='x (Angstroms)',
                    yaxis_title='y (Angstroms)',
                    zaxis_title='z (Angstroms)'),
                    width=700,
                    margin=dict(r=10, b=10, l=10, t=10))

fig.show()
In [ ]: